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Abstract. This work concerns the development of an Algebraic Multilevel method 
for computing stationary vectors of Markov chains. We present an efficient Boot- 
strap Algebraic Multilevel method for this task. In our proposed approach, we 
employ a multilevel eigensolver, with interpolation built using ideas based on com- 
patible relaxation, algebraic distances, and least squares fitting of test vectors. Our 
adaptive variational strategy for computation of the state vector of a given Markov 
chain is then a combination of this multilevel eigensolver and associated multilevel 
preconditioned GMRES iterations. We show that the Bootstrap AMG eigensolver 
by itself can efficiently compute accurate approximations to the state vector. An 
additional benefit of the Bootstrap approach is that it yields an accurate inter- 
polation operator for many other eigenmodes. This in turn allows for the use of 
the resulting AMG hierarchy to accelerate the MLE steps using standard multigrid 
correction steps. Further, we mention that our method, unlike other existing mul- 
tilevel methods for Markov Chains, does not employ any special processing of the 
coarse-level systems to ensure that stochastic properties of the fine-level system are 
maintained there. The proposed approach is applied to a range of test problems, 
involving non-symmetric stochastic M-matrices, showing promising results for all 
problems considered. 



1. Introduction 
We consider the task of computing a non-zero vector x such that 
(1) Ax = 1 X, 

where A denotes the transition matrix of a given irreducible Markov process and x is 
the associated state vector, an eigenvector of A with eigenvalue equal to one. Since 
A is irreducible, x is unique up to scalar factors. The approach considered in this 
paper is to approximate x iteratively using an Algebraic Multigrid (AMG) method 
designed to find a non-trivial solution to the eigenvalue problem 

Bx = X, where B = I — A^ 

in the setup phase and use the computed AMG method to solve the homogeneous 
system 

Bx = 0, 

. Our AMG algorithm relies on the Bootstrap framework, a fully adaptive algebraic 
multigrid scheme proposed in [4j. For solving complex- valued linear systems, this 
Boot AMG (Bootstrap AMG) framework was efficiently put into action in [6]. In this 
paper, we develop a variant of Boot AMG specifically tailored to compute x in ([T]). 
Although we are not (yet) in a position to give a full rigorous mathematical analysis 
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of the method, we demonstrate the efficiency of the BootAMG approach for a series 
of Markov Chain test problems. 

The main new ingredient of our BootAMG method is an adaptive (setup) compo- 
nent that simultaneously computes an approximation to the state vector and Test 
Vectors (TVs) which are then used in a Least Squares approach to calculate a multi- 
level sequence of interpolation operators. Letting ^ denote a given level of the multi- 
level hierarchy and taking on the finest level Bq = our adaptive strategy consists 
of the following steps (explained and described in detail in the latter sections): 

(1) relax on the homogeneous system B^x^ = 0,S^ G M^^^^^ to improve the set 
of test vectors, xf\...^xf\ approximating the state vector; rii here is the 
problem size on level £; 

(2) compute a sparse interpolation operator G W^w^^rn using a least-squares 
fit based on x^l\ ...^x^p] 

(3) compute the system matrix S^+i G ]R^^+ix^^+i on the next coarser level as 
BiJ^i = QiB^Pi via a Petrov-Galerkin approach, with the restriction Qi G 
jgn^xn^+i representing a simple averaging operator; also calculate the mass 
matrix T/+i = QiTiPi (with Tq = /) for the next leg of BootAMG; 

(4) obtain the coarse-level series of test vectors 

= R£x[^\ j = 1, with Ri G being injection. 

This process is applied recursively until a level with sufficiently small problem size 
(the coarsest level, ^ = L) is reached. There, an exact eigensolve for 

BiX — XT^x 

is performed and then used to improve the given approximations to the state vector 
on increasingly finer levels. Along with the state vector's approximations x^ ^ some 

(small number) of near kernel eigenfunctions (A^,x^^^), from the coarsest level are 
interpolated to increasingly finer levels, where they are relaxed and then used as 
initial guesses to the corresponding fine-level eigenvectors. On each such level, the 
relaxation is performed on the homogeneous system 

{B, - \iT,)xf = 0, 

and the approximations to are updated (except for Aq which is known to be exactly 
zero.) Each x^ is already a reasonable approximation to an eigenfunction on level £, 
i.e, 

B,xf ^ X%xf 

for an i^^ level eigenvalue Af. Therefore, a good eigenvalue approximation can be 
obtained by 

{TiX^.v) 

for any vector v. We choose v — x^ and thus calculate an approximation to eigen- 
value using what would be a Rayleigh quotient for Hermitian matrices. The change in 
eigenvalue approximation for each vector x^'^\ i > 1 after relaxation provides a mea- 
sure of their accuracy - if after relaxation the relative change is significant, we add 
x^^^ to the set of test vectors. Otherwise, these eigenvectors are well approximated by 
the existing interpolation operator to the current level and thus need not be built into 
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the related coarser space. The state eigenvector approximation is always included in 
the test set. After all test sets on all levels, including the finest, are updated, all P^, 
and are recalculated, and the setup stage is complete. 
Once the adaptive setup is finished, we use the resulting V-cycle AMG method 
as a preconditioner to (full) GMRES, applied to the homogeneous system Bx = 0. 
This iteration converges rapidly to the state vector, since as we demonstrate later, the 
Boot AMG preconditioner efficiently separates the state vector from the field of values 
of B on an appropriate complementary subspace. This use of both MLE steps and 
V-cycle preconditioned GMRES iterations is the other main new idea in the present 
paper. 

The remainder of this paper is organized as follows. Section [2] contains an introduc- 
tion to Markov chain systems and a general review of AMG approaches for computing 
state vectors of Markov chains. Section |3] contains a general description of the Boot- 
strap AMG ideas along with full algorithmic descriptions of realizations that we use 
in this paper. In Section [4] we discuss our BootAMG algorithm for Markov chains 
and theory for our AMG-GMRES solver. Numerical experiments are presented in 
Section [5] and concluding remarks are given in Section [6| 



The transition matrix A G W^^^ of a Markov process contains as its entries a^j, the 
transition probabilites from state i to state j, a^j > for all i^j. Matrix A is column 
stochastic, i.e. A^l = 1, with 1 being the vector of all ones. It is always possible to 
eliminate self-transitions [19j, so we assume an = for all i from now on. The state 
vector X satisfies 



with X 7^ 0, < x^. By the Perron-Frobenius theorem (cf. [1^ p. 27], e.g.) such a 
vector X always exists. 

A general square matrix A G M.^^^ induces a directed graph G{A) = (y, E) with 
vertices 1/ = {1, . . . ,n} and directed edges E = {{i^j) C : i j and a^j 0}. 
Two vertices i and j in G{A) are said to be strongly connected (as in graph theory) 
if there exist directed paths in G{A) from i to j and from j to i. Since this is an 
equivalence relation on the vertices, it induces a partitioning of V into the strong 
components of G{A). If G{A) has exactly one strong component, the matrix A is 
called irreducible; otherwise it is called reducible which is precisely the case when 
there exists a permutation matrix V such that 



where An and A22 are square submat rices. 

If A is irreducible — which we assume throughout — the Perron-Frobenius Theorem 
guarantees that the state vector x has all positive values and is unique up to a 
scalar factor. From this theorem it also follows that the spectral radius of A satisfies 
^(^4) = 1, implying that B = I — A is Sl singular M-matrix; recall that < a^j < 1. 
In addition, B is irreducible, because of the irreducibility of A. 

The idea of using Multigrid (MG) to compute the state vector of an irreducible 
transition matrix is not new. Numerous approaches have been explored in the past, 
in both the aggregation-based [lOlIIIlEOlEI] and Classical [23] MG settings. In [20], 
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the idea of using a partition of unity type aggregation method was first explored. 
Later, in [TOj, an extension of this approach based on Smoothed Aggregation was 
developed for application to the page rank problem and further generalized in pTj. 
For an overview of Classical MG type methods (and preconditioners) we direct the 
reader to [23] and the references therein. 

Our MG solver can be loosely categorized as a variation of the Classical AMG 
approach in the following way. We choose the coarse level variables as a subset of the 
fine level variables. In this aspect, our approach is most closely related to the recent 
work by Virnik [23j; we too employ the AMG-type method as a preconditioner for 
GMRES. Our proposed new approach, however, differs from all previous AMG-type 
solvers for Markov chains in several ways. Most importantly, we build interpolation 
adaptively using a Least Squares approach, and we combine the multilevel eigensolver 
with application of correction steps using an AMG-GMRES preconditioner. Unlike 
other methods, we thus (experimentally) observe an excellent convergence without 
the need to recompute interpolation and the entire coarse-level hierarchy at each 
iteration. In addition, our solver yields an efficient and straightforward method for 
Markov chains which computes the state vector to arbitrary accuracy without the 
need for lumping (see [llj) and other such approaches typically employed to maintain 
the stochastic properties of the transition matrix and state vectors on all coarse levels 
of the MG hierarchy. 

3. Bootstrap AMG 

In this section, we provide a general review of the Bootstrap AMG framework [HIT] 
for building MG interpolation together with some heuristic motivation of our choices 
for the individual components of the BootAMG multigrid algorithm. 

The first key ingredient to any AMG method is its relaxation or smoothing iteration. 
For a given system matrix 5 it is a splitting based iteration 

(2) Mu""^^ = Nu"" + 6, where B = M - N, M non-singular, = 0, 1, . . . 

with its "error propagation matrix" given by M~^N = I — M~^B. Here, M is 
chosen such that after a few iterations the error w.r.t. the solution of Ex = 6 is 
algebraically smooth^ i.e. ||5e^|| <C In many situations this can be achieved by 

point-wise Gauss-Seidel or (under-relaxed) Jacobi iterations. Indeed, in our Markov 
chain setting we used cj- Jacobi relaxation with uj = 0.7 on all levels, i.e., we took 
M = ^diag(S), resulting in / - M'^B = I - u - {didig{B))-^ B . 

(In case a diagonal entry ba of B happens to be zero for some row(s) i, Kaczmarz 
or some other distributive relaxation scheme can be applied to the i^^ equation.) 

The Bootstrap AMG process can now be described as follows. Coarse variables 
are selected as a subset of the fine variables using a full coarsening (for problems on 
structured grids) or, otherwise, compatible relaxation (CR) coarsening scheme [3j. 
The CR scheme can be either started from scratch, or, if geometric information is 
given and a suitable candidate set of coarse variables is known, such set can be tested 
and improved by CR. Interpolation is then computed using a Least Squares based 
approach. We mention that once a tentative Multigrid hierarchy has been defined, it 
can be used to further enhance the set of test vectors. 
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3.1. Choosing the coarse variables: compatible relaxation. In the AMG set- 
ting, CR is a relaxation-based coarsening process which can be viewed as a special 
case of a general approach for constructing coarse-level descriptions of given fine-level 
systems, including non-stationary, highly non-linear, and also non-deterministic sys- 
tems [5j. The basic idea of CR is to use the given relaxation scheme ([2]), restricted to 
appropriately defined subspaces, to measure the quality of the given coarse space and 
also to iteratively improve it if needed. We proceed with a brief overview of CR and 
its use in AMC coarsening. A detailed discussion, theory and comparisons between 
various measures of the quality of coarse spaces and their relations to compatible 
relaxation schemes are presented in [3l El O [121 1131 [IS] . 

3.1.1. Classical AMG CR-based coarsening. Assume that the set of coarse-level vari- 
ables, C, is a subset of the set of fine-level variables, Under this assumption, one 
possible form of CR is given by J^-relaxation for the homogeneous system — relax- 
ation applied only to the set of variables, with := Q\C. Civen the partitioning 
of Q into J-" and C, we have 

" - [uj ' ^ - [B^f B^J ' M - ^^^^ ^^J , 

assuming the equations are permuted such that the unknowns in come before those 
in C. The J^-relaxation of CR is then defined by 

(3) u^+' = {I-Mf}BffX^ = EfU^f. 
If M is symmetric, the asymptotic convergence rate of CR 

where p denotes the spectral radius, provides a measure of the quality of the coarse 
space, that is, a measure of the ability of the set of coarse variables to represent error 
not eliminated by the given fine-level relaxation. This measure can be approximated 
using J^-relaxation for the homogeneous system with a random initial guess t^j. Since 

lim^^oo 11^/11''^^^ = P{^f) iiorm || • ||, the measure 

(4) {\\u)\\m\\Y'^ 

estimates pf. 

In choosing C, we use the CR-based coarsening algorithm developed in [9j. This 
approach is described in Algorithm [TJ 

In our numerical experiments for Markov chains, we use weighted Jacobi CR, set 
the CR tolerance = .85, the number of CR sweeps = 8, choose the components of 
to be uniformly distributed in the interval [1,2], and select Co using the standard 
independent set algorithm (see [22j) based on the full graph of the system matrix. 
For further information on CR we refer the reader to [9j . 

3.2. Building Bootstrap AMG interpolation. We now outline the Least Squares 
(LS) approach for defining interpolation, see also [6j. We assume that the set C of 
coarse variables is given, i.e., it has been previously determined, for instance, by 
geometric coarsening and/or using compatible relaxation. In the LS approach, the 
interpolation operator P is chosen to approximate a given (specifically chosen) set 
of test vectors. We define c as the maximum number of coarse-level variables used 
to interpolate to a single fine-level variable, or equivalent ly the maximum number of 
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Algorithm 1 compatible_relaxation {Computes C using Compatible Relaxation} 

Input: Co {Co = is allowed} 
Output: C 
Initialize C = Co 
Initialize J\f = ^}\C 

Perform u CR iterations ^ with components of randomly generated 
while pf > 9 do 

Af = {ien\C:j^^>9} 

I'^i I 

C = C U {independent set of JV} 

Perform u CR iterations ^ with components of randomly generated 
end while 



non-zero entries in any row of P. The key ingredient of the BootAMG setup lies in 
the use of several test vectors (TVs) that collectively should represent those error 
components not reduced by relaxation. We assume for now that such a set of test 
vectors, U := {x^^^l^^ is known on some fine level. The rows of the prolongation 
operator P are then obtained individually. For each variable i G J-", we first determine 
a set of its neighboring coarse-level variables, A/^, using the directed graph G{A) 

(5) A/i = {j G C : there is a path of length < z in G{A) from i to j }, 

where z < 3, such that Afi is a subset of the local graph neighborhood of i. We 
then determine an appropriate set of (coarse level) interpolatory variables C J\f^ 
with \J^i\ < c that yields the best fit of the LS interpolation for point i ^ T . For 
ease of presentation let us use the indices j G Ji of the fine level variables to address 
the columns of P . We then define the local least squares functional for the nonzero 
entries Pi = {pij : i G JT^} of row i of P as 

(6) m; = j2^k U'^ - E Prjx^ 

k=o \ jeJi 

The task is then to find a set J^i of interpolating points for which the minimum of 
L is small and to obtain the corresponding values Pij of the minimizer that yield the 
coefficients for the interpolation operator. Generally, the weights should be chosen 
according to the algebraic smoothness of x*^^^ to bias the least squares functional 
towards the smoothest vectors. We do so by using the square of the inverse of the 
norm of the residual of the homogeneous problem B^x^^^ = for each x*^^^.Here £ is 
the current level to which operator P interpolates from the next coarser level £ + 1. 

Our approach for solving this task is given as Algorithm [2} It is based on a greedy 
strategy to find an appropriate set of interpolatory points J^i for each fine point i] see 
[6] and [IS]. 

4. MLE AND AMG-GMRES Algorithms 

In this section, we discuss the main ideas and their implementations in our Boot AMG- 
based algorithm for computing non-trivial solutions x to Bx = 0, where as before 
S = / — ^4, ^4 is an irreducible column stochastic transition matrix so that S is a 
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Algorithm 2 ls_interpolation {Computes Least Squares based interpolation} 

Input: c 

Output: interpolation P 

for i G do 

Take A/^ from g 
Set J, = 
repeat 

determine g"" G JVi s.t. minL(P^; Ji U {s'*}) = minminL(P^, J"^ U {g}) 

Pi geAfi Pi 

Set A/; = A/; \ {/} and J, = J, U {/} 
until IJ^I > c or A/'i = 
end for 



non-symmetric and singular M-matrix. The goal of each MLE step is twofold: (1) 
to compute approximations to the solution x and also to compute suitable test vec- 
tors that allow us to improve the accuracy of Least Squares interpolation and (2) to 
update the multigrid hierarchy, used to precondition GMRES. 

Such a hierarchy consists of a sequence of prolongation and restriction operators 
and systems of coarse-level equations. The coarse-level equations are defined using 
the Petrov-Galerkin approach. More precisely, S^+i = QiBiPi^ with restriction Qi 
given by the averaging operator. This means that each column of Qi has identical 
nonzero entries at exactly the positions corresponding to Ji and that l^Qi = 1^ This 
implies Bjl = for all levels /. 

4.L Multilevel Eigensolver (MLE). Our approach for building a multilevel method 
relies on the BootAMG framework. We construct a sequence of increasingly coarser- 
level descriptions of the given fine-level Markov Chain system using the LS-based 
fit discussed above. The coarse-level test vectors are calculated from the fine-level 
test set as X^^^ = Rx^^\ On the coarsest level, the eigenproblem is solved exactly 
for some small number of the lowest (in absolute value) eigenvectors, which are then 
directly interpolated and relaxed on finer levels, becoming increasingly more accurate 
approximations to the finest-level eigenfunctions (such treatment is similar to the Ex- 
act Interpolation Scheme, e.g, fTT]). These eigenfunction approximations are added 
to test sets on each level. Note that the initial (random) test functions are not relaxed 
during this stage - otherwise they may become almost linearly dependent which may 
result in ill-posed local least squares problem. With the improved TV sets (enriched 
by eigenfunction approximations), new prolongation operators are then calculated. 

We note that the setup cycle in Algorithm [3] also yields an approximation to the 
state vector which we use as an initial guess in our AMG-GMRES iterations. This 
approach can also be used as a stand-alone solver for these systems although, in 
general, it will be far more expensive than using several AMG-GMRES steps in 
addition to the MLE. We note that for certain problems (e.g., when several eigenvalues 
of B are close to zero) and for larger problem sizes, it may be necessary to alternate 
several times between the MLE and AMG-GMRES iterations to find an accurate 
approximation of x. In our tests, we found it sufficient to use either a traditional 
l/-cycle (/i = 1) or Vl^-cycle (/x = 2) for the MLE step followed by a small number of 
AMG-GMRES steps to compute an approximation to the state vector. Ultimately, 
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in our MLE approach, the sequence of prolongation operators needs to be accurate 
for only the smoothest error component - the kernel of the finest-level operator B = 
I— A. The accurate resolution of many other near- kernel components by our LS-based 
interpolation operator allows for fast convergence of MG preconditioned GMRES, 
whereas the kernel component of B is the state vector we aim to compute. Algorithm 
[3] contains a pseudo-code of our implementation of the MLE iteration. 



Algorithm 3 bootamg_mle {Implements the Bootstrap AMG MLE scheme} 

Input: Bi {Bq = B)^ Ti^ (Tq = I)^ Ui (set of test vectors on level /) 
Output: updated Ui and (A/, V^), approximations to the lowest eigenpairs; 
if / = L then 

Compute Vl = {x^^}i,...,k^ s.t. S^x^ = A^r^xW, \X^^^\ < \X^^^\ < . . . < |A(^)|. 
else 

Relax Bix^^^ = 0, x^^^ G Ui 
Set V/ = 

for m = 1, . . . , /i do 

Pi = ls_interpolation(Z/// U V/, c) 

Compute averaging Qi with sparsity((50 = sparsity(P/) 

Bi^i = QiBiPi 

Ti^i = QiTiPi 

Ui+i = {Rix.x G Ui} 

Vi+i = bootamg_mle(5/+i,r/+i,Z//^+i) 

Vi = {Pix,xeVi^i} 

for i = 1, . . . , |V/| do 

Relax {Bi - \iTi)x^^ = 0, x» G Vi 

end for 
end for 
end if 



In Figure [T] a possible setup cycle is visualized. We note that at each square one has 
to decide whether to recompute P or advance to the next finer grid in the multigrid 
hierarchy. The illustrated cycle resembles a Vl^-cycle, but any other cycling strategy 
can be applied to the setup. 

4.2. Multigrid preconditioned GMRES. The multilevel eigensolve steps described 
above yield increasingly better approximations to the state vector and other vectors 
that cannot be removed efficiently by the Multigrid relaxation. 

Although we are only interested in computing the state vector, the LS-based Multi- 
grid hierarchy is able to resolve a larger subspace. This leads to the idea of exploiting 
this richness of the given hierarchy for use in MG correction steps - in addition to 
the discussed MLE steps. 

To illustrate the effect of MG correction steps applied to the homogeneous prob- 
lem Bx = 0, we start by analyzing simple relaxation schemes for the steady state 
problem, Ax = x. As the steady state solution is the eigenvector corresponding to 
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• Relax on Bu = 0,u eU • Compute V, s.t., Bv = XTv, v eV 

• Relax on {B - \T)v = 0, G V Relax on Bu = {),u eU and {B - \T)v = 0, G V 

Figure 1. Bootstrap AMG setup W-cycle. 

the eigenvalue with largest absolute magnitude, a power iteration 

(7) = Ax^ 

is guaranteed to converge to the solution. However, convergence can be slow if A has 
other eigenvalues with absolut value close to one. 

Such power iterations applied to the steady state problem are in turn equivalent to 
applying a Richardson iteration to the homogeneous system [I — A)x = Q. A natural 
modification, facilitating that the field of values of A is contained in the unit circle, 
is then given by a suitable under-relaxed iteration, yielding the error propagator 

(8) e^+^ = {I-TB)e\, 
which translates into an identical relation for the iterates, 

(9) = {I-r{I-A))x\ 

so that the iteration can be interpreted as a modified power method. 
In Figures 2(a) and 2(b)| the spectra of A and of / — r (/ — ^4) are depicted for a 



characteristic two-dimensional test problem, along with the field of values. The error 
propagator of the MG V-cycle applied to Bx = is given by 

(10) E = {I - MB) (/ - PE^QB) (/ - MB) , 

where denotes the error propagation operator on the next coarser level. This 
recursive definition terminates at the coarsest level with E = B\ the Moore- Penrose 
inverse of the matrix on the coarsest level, since we compute the minimal norm 
solution of the respective system at that level. The error propagator can be rewritten 
as 

(11) I-CB, 

so that comparing with the one from ([s]) (relaxed Richardson) we see that we can 
interprete the multigrid V-cycle as yet another modification of the power iteration. In 



Figure [2(c) , the spectrum and field of values of the MG V(2, 2)-cycle preconditioned 
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(a) Spectrum and field of val- (b) Spectrum and field of val- (c) Spectrum and field of val- 
ues for Richardson error prop- ues for r-Richardson. ues of Multigrid error propa- 
agator. gator. 



Figure 2. Spectra and field of values of the Richardson, r-Richardson 
and Multigrid error propagators for two-dimensional tandom queuing 
problem with n = 33^. 



matrix (11) is depicted. For this case, it is clear that applying power iterations to 
this preconditioned operator will converge rapidly to the steady state vector. 

To further accelerate this approach, instead of the straightforward power method, 
we consider MG preconditioned GMRES steps. As we show below, such iteration is 
guaranteed to converge to the solution of our Markov Chain systems. 

Remark 4.1. Let A be a column- stochastic irreducable operator and B = I — A. 
Then we have 

(12) range {B) n null{B) = {0}. 

Proof. As A is column-stochastic we know that A^l = 1. Hence B^l = 0. Fur- 
thermore we know from the Perron-Frobenius theorem that Bx = is uniquely 
solvable up to a scalar factor and the solution x* is strictly positive (component- 
wise). Thus, as for all y G range (5) there exists a z such that y = Bz we have 
(1,2/) = {l,Bz) = {B^l.z) = 0. With this and (l,x*) as X IS strictly positive 



we get (12). □ 



In [m Theorem 2.8] it is shown that GMRES determines a solution of Bx 



for all b G range (5) ,Xo G iff range (S) n null (5) = {0}. Due to Lemma 4.1 
the assumptions of this theorem are fulfilled for B = I — A^ where A is a. column- 
stochastic irreducible operator arising in Markov-Chain models. There seems to be 



no simple way to extend Lemma |4.1| to the multigrid preconditioned matrix. So the 
observation that GMRES for the multigrid preconditioned matrix worked out very 
well in all our experiments is not backed by theory so far. 

4.3. A multigrid preconditioned Arnoldi method. Instead of solving the homo- 
geneous linear system with preconditioned GMRES, we can also perform the Arnoldi 
method for the preconditioned matrix CB. We thus compute orthonormal vectors 
Vi^V2^ . . . in the usual way such that Vi^ . . . is a basis of the Krylov subspace 
span{^, C?/, . . . , (CS)^~-^y}. Denoting G M^^^ the orthogonal projection of CB 
onto the Krylov subspace, i.e. = V^{CB)Vk with Vk = . . . 1'^^], we compute 
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the eigenvalue closest to of and the corresponding eigenvector rj. We then take 
the Ritz vector VkTj as an approximation for the steady state vector. For the vector y 
upon which we build the Kryolv subspace we take the approximation x for the steady 
state vector resulting from the MLE setup phase. 

5. Numerical Results 

In this section, we provide results obtained using our Bootstrap AMG method when 
applied to a series of Markov Chains. Our numerical tests consist of our BootAMG- 
based approach to three Markov chain models that can be represented by planar 
graphs. Each of the models has certain interesting characteristics that pose problems 
for the solver. 

5.1. Test Problems. We begin our experiments with a very simple model. The 
uniform two-dimensional network can be seen as the Markov chain analogue of the 
Laplace operator. It is defined on an equidistant grid Vt of size N x N. We denote 
this grid in graph notation as Gq = (Vq^ Eq). The entries of A are then given as 

^'■^ \ 0, otherwise, 

where dout{j) is the number of outgoing edges oi j ^ Q. In Figure |3| we illustrate the 
two-dimensional uniform network problem. In the tests we conduct, we use again full- 
coarsening, i.e., we choose C to be every other grid point in both spatial dimensions. 
In order to keep the overall invested work small, we consider to take only up to two 
interpolatory points Ci per point i G = Q\C. As the steady-state vector is known 
to be strictly positive, we choose to use initially random, but positive test vectors for 
this first problem, and also in all following tests. More precisely, we choose vectors 
with entries uniformly distributed in [1,2]. 

In Table [!} we present results that use 6 test vectors and a V{2, 2)-cycle MLE 
step with cj-Jacobi smoother, u) = .7. We report the number of iterations needed to 
compute the steady-state vector x, such that 

\\Bx\\ < 10"^ ||x|| = L 




Figure 3. Uniform network model on a two-dimensional equidistant grid. 




Figure 4. Tandem-Queuing Network with probability to advance /x 
and queuing probabilities /Hx and /Hy on a two-dimensional equidistant 
grid. 

In addition, the number of preconditioned GMRES iterations needed to achieve the 
same accuracy is reported, with the initial MLE setup cycle denoted by the subscript. 
The coarsest grid in the experiments is always 5x5. The operator complexity in these 



N 


17 


33 


65 


129 


MLE 


10 


9 


10 


11 


pGMRES 


7v 


Sv 


lOy 


10y2 


pArnoldi 


7v 


Sv 


lOy 


9l/2 



Table 1. Multilevel results for the two-dimensional uniform network 
model on an X grid. We report results to compute the steady- 
state vector X to an accuracy of 10~^, using a V{2^ 2)-MLE cycle with 
cj-Jacobi smoothing, uj = .7. In addition we also report the number 
of iterations pGMRES needs to achieve the same accuracy, where we 
denote the initial bootstrap setup in the subscript. The sets U and V 
consist of 6 initially positive random vectors and coarsest grid eigen- 
vectors, respectively. 



tests is bounded by 1.6. The test shows that the MLE scales with increasing problem- 
size, while the number of preconditioned GMRES iterations grows slightly. However, 
one step of pGMRES is much cheaper than one step of MLE. Although we do not 
report results here, in general we may also consider restarting the MLE iterations 
when the solution is not found in a reasonable number of pGMRES iterations. 

The next Markov chain model we consider in our tests is a tandem queuing network, 
illustrated in Figure |4} The spectrum of this operator is complex, as we show in 
Figure |2(a)[ Again, we use full-coarsening and a coarsest grid of 5 x 5. We present 
our results in Table [2} In the test we use the following probabilities, 

_ 11 _ 10 _ 10 

^"31' ^""31' ^'~31' 
Again, we see that the MLE method converges rapidly to the steady-state vector 
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17 


33 


65 


129 


MLE 


8 


8 


8 


8 


pGMRES 


6y 


6a/ 


6y 


7y 


pArnoldi 


6y 


6y 


6y 


7y 



Table 2 . Multilevel results for the tandem queuing network model on 
an X grid. We report results to compute the steady-state vector 
X to an accuracy of 10~^, using a ]/(2,2)-MLE cycle with cj-Jacobi 
smoothing, u = .7. In addition we also report the number of iterations 
pGMRES needs to achieve the same accuracy, where we denote the 
initial bootstrap setup in the subscript. The sets U and V consist 
of 6 initially positive random vectors and coarsest grid eigenvectors, 
respectively. 




Figure 5. Approximations to the eigenvectors corresponding to the 
smallest 6 eigenvalues of the tandem queuing network problem on a 
129 X 129 grid with 6 levels upon convergence of the steady-state solu- 
tion of the MLE method. 

and also yields a very efficient preconditioner for the GMRES method. We observe 
that the number of iterations does not depend on the size of the problem. Note 
that the MLE method also yields accurate approximations to the eigenvectors corre- 
sponding to all k smallest eigenvalues on the finest grid. In Figure [5| we show the 
computed approximations and report in Table [3] their accuracy upon convergence of 
the steady-state vector. One should keep in mind that the results are intended 



i 


1 


2 


3 


4 


5 


6 




1.02£;-8 


9.04£'-3 


2.68£;-2 


2.84£'-2 


1A2E-1 


3.77£'-l 



Table 3. Accuracy of the eigenvectors corresponding to the smallest 
6 eigenvalues of the tandem queuing network problem on a 129 x 129 
grid with 6 levels upon convergence of the steady-state solution of the 
MLE method. 



to show the promise of our approach rather than presenting an optimized method 
in the bootstrap framework for this type of problems. We limit our analysis to the 
statement that with minimal effort spent in adjusting the parameters of the setup of 
the bootstrap approach we obtain scalable solvers. The optimization of the method 
by tuning the parameters involved, e.g., relaxation parameter of the smoother, coars- 
ening, caliber, number of test vectors, weighting in the least squares interpolation, 
number of relaxations in the setup and solution process, is part of future research. 
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The last test we consider corresponds to a triangulation of randomly chosen 
points in [0, 1] x [0, 1]. The transition probabilities in the network are then given by 
the inverse of the number of outgoing edges at each point, similar to the uniform 
network. In Figure [6} we show two examples of such networks, one with = 256 
and one with = 1024 points. Due to the fact that the corresponding graphs of this 




Figure 6. Delaunay-triangulations of randomly chosen points in 
the unit square [0, 1] x [0, 1]. 

model are planar, we call this model unstructured planar graph model. As there is no 
natural way to define the set of coarse variables C for this problem we use compatible 



relaxation, introduced in section [3T] to define the splitting of variables Q = U C. 
In Figure [7| a resulting coarsening is presented. The size of each individual point 
represents how many grids it appears on. In Table [4| we report results of our overall 
algorithm for unstructured planar graph models for a set of varying graphs. Even 



N 


256 


512 


1024 


2048 


MLE 


15 


20 


20 


20 


pGMRES 


Sv 


lOy 


lOy 


lly 


pArnoldi 


Sv 


lOy 


lOy 





TABLE 4. Multilevel (3 level) resu! 



ts for the unstructured planar graph 



model on a grid with A^ points randomly scattered in the unit square. 
We report results to compute the steady-state vector x to an accuracy 
of 10~^, using a 1/(2, 2)-MLE cycle with cj-Jacobi smoothing, uj = .7. 
In addition we also report the number of iterations pGMRES needs to 
achieve the same accuracy, where we denote the initial bootstrap setup 
in the subscript. The sets U and V consist of 6 initially positive random 
vectors and coarsest grid eigenvectors, respectively. 



for this unstructured graph network we obtain a fast converging method with our 
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(a) Random planar graph with N = 256 (b) Random planar graph with N = 1024 



Figure 7. Coarsening of Delaunay-triangulations of points ran- 
domly scattered in the unit square [0, 1] x [0, 1] shown in Figure |6] using 
compatible relaxation. 

MLE approach. The mild variation in the results when increasing the problem-size 
might be caused by the fact that by increasing the problem-size the nature of the 
problem changes in the sense that the average number of outgoing edges of grid 
points increases. That is, it is not necessarily clear whether two unstructured graphs 
of different sizes are comparable. 

6. Conclusions 

The proposed Boot AMG Multilevel setup algorithm produces an effective multilevel 
eigensolver for the Markov-Chain test problems we considered. In general, we have 
proposed and implemented two important new ideas: First, the use of a coarse-level 
eigensolve and the resulting multilevel hierarchy to improve a given approximation of 
the state vector. Second the use of Boot AMG preconditioned GMRES steps to further 
accelerate this eigensolver. Both ideas, separately or combined can be incorporated 
into any given multilevel method used for solving Markov Chain problems or other 
problems targeting smooth eigenvectors. The accurate representation of the near- 
kernel of the fine-level system on coarser levels, that results from our BootAMG setup, 
yields a very effective preconditioner to GMRES as well as for the Arnold! method. 
An additional benefit of our proposed method over other existing multilevel method 
for Markov Chains is that we do not require any special processing of the coarse-level 
systems to ensure that stochastic properties of the fine-level system are maintained 
there. We mention that the developed approach is not restricted to Markov Chain 
problems; it can be applied to other eigenvalue problems. 
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